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Abstract — In this paper, a novel maximum Doppler spread 
estimation algorittim is presented for OFDM systems with the 
comb-type pilot pattern in doubly selective fading channels. First, 
the least-squared estimated channel frequency responses on pilot 
tones are used to generate two auto-correlation matrices with 
different lags. Then, according to these two matrices, a Doppler 
dependent parameter is measured. Based on a time-varying 
multipath channel model, the parameter is expanded and then 
transformed into a non-linear high-order polynomial equation, 
from which the maximum Doppler spread is readily solved by 
using the Newton's method. The delay-subspace is utilized to 
reduce the noise that biases the estimator. Besides, the subspace 
tracking algorithm is adopted as well to automatically update the 
delay-subspace. Simulation results demonstrate that the proposed 
algorithm converges for a wide range of SNR's and Doppler's. 

Index Terms — Doppler spread. Estimation, Subspace tracking, 
OFDM, Doubly selective fading channels. Comb-type pilot. 

I. Introduction 

The maximum Doppler spread is one of the key parameters 
for the adaptive strategies to tune mobile communication 
systems to accommodate various radio transmission environ- 
ments [1] and, especially, alleviate the inter-carrier interfer- 
ence (ICI) for orthogonal frequency division multiplexing 
(OFDM) systems. Most methods of estimating the maximum 
Doppler spread reported in literatures are categorized into two 
classes [2]: level-crossing-rate-based and covariance-based 
techniques. Since the algorithms reviewed in [2], such as [3], 
were not specifically designed for OFDM systems, they did 
not exploit the special signal structure of OFDM systems. 
On the other hand, most algorithms designed for OFDM 
systems are covariance-based. [4] determined the maximum 
Doppler spread through estimating the smallest positive zero 
crossing point. Cai [5] proposed to obtain the time correlation 
function (TCF) by exploiting the cyclic prefix (CP) and its 
counterpart. However, Yucek [6] commented that for scalable 
OFDM systems whose CP sizes were varying along the time, it 
was difficult for [5] to sufficiently estimate the TCF, therefore 
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its accuracy would be degraded significantly. In stead, Yucek 
proposed to estimate the maximum Doppler spread by taking 
advantage of the periodic training symbols. However, in order 
to reduce overheads, training symbols are sparse and typically 
transmitted as preambles to facilitate the frame timing and 
carrier frequency synchronization, which would cause [6] to 
converge slowly or even fail. 

In this paper, we propose to estimate the maximum Doppler 
spread by exploiting the comb-type pilot tones [7] that are 
widely adopted by wireless standards. The channel frequency 
responses (CFR's) estimated from the pilot tones are projected 
onto the delay-subspace [8] to reduce the noise perturbation 
and then used to acquire a Doppler dependent parameter. With 
a careful expansion of the parameter, a nonlinear high-order 
polynomial equation is formed, from which the maximum 
Doppler spread is readily solved by resorting to the Newton's 
iteration. Moreover, the subspace tracking algorithm [9] is 
adopted as well to track the drifting delay-subspace. 

This paper is organized as follows. In Section HH the OFDM 
system and channel model are introduced. Then, the maximum 
Doppler spread estimation algorithm is presented in Section 
Uni Simulation results and analyses are provided in Section 
HVl Finally, Section IV] concludes the paper 

A. Basic Notation 

Uppercase and lowercase boldface letters denote matrices 
and column vectors, respectively. (•)^ and || • \\p denote 
conjugate transposition and Frobenius norm, respectively. E{-) 
represents the mathematical expectation of a stochastic pro- 
cess, and denote the i-th and (i,j)-th elements of a 
vector and a matrix, respectively. diag{A) denotes a diagonal 
matrix with the diagonal elements of A on the diagonal. 

II. System Model 

Consider an OFDM system with a bandwidth of BW = 
l/T Hz (T is the sampling period). N denotes the total number 
of tones, and a CP of length Lcp is inserted before each 
symbol to eliminate the inter-block interference. Thus, the 
whole symbol duration is Tg ~ {l+rcp)NT, where r, '""^ 



' cp 



N 



In each OFDM symbol, P tones are transmitted as pilots to 
assist channel estimation. In addition, the optimal pilot pattern, 
i.e., the equipowered and equispaced [10], is assumed. 

The complex baseband model of a linear time-variant mo- 
bile channel with L paths can be described by [11] 



h{t,T) 



^1=0 



hi{t)5 {t - TiT) 



(1) 



2 



where ti E TZ is the non-sample-spaced delay of the ^-th path 
normalized to the sampling period, and hi{t) is the corre- 
sponding complex amplitude. According to the assumption 
of the wide-sense stationary uncorrected scattering, /i;(t)'s 
are modeled as uncorrected narrowband complex Gaussian 
processes. In the sequel, P > Lis assumed for determinability. 
Furthermore, by assuming the uniform scattering environment 
introduced by Clarke [12], all /i;(t)'s have the identical nor- 
malized TCF, therefore the TCP of the l-th path is 



PxL 



as 



rt,iiAt) =af Jo {27: fdAt) 



(2) 



where af is the power of the l-th path, fd is the maximum 
Doppler spread, and Jo( ) is the zeroth order Bessel function 
of the first kind. Moreover, the total power of the channel is 
normalized, i.e., '^^Sq (jf — 1. 

Assuming a sufficient CP, i.e., L^p > L, the discrete signal 
model in the frequency domain is 



y/(n) = H/(n)x/(n) + n/(n) 



(3) 



where x/(n), y/(n), n/(n) S qNxi j-j^g ^ transmit- 
ted and received signal and additive white Gaussian noise 
(AWGN) vectors, respectively, and B.f{n) e Q^^^ is the 
channel transfer matrix whose (fc + i', fc)-th element is 



[H/(n)] 



k-\-v,k 



N-lL-1 
m=0 1=0 



n, mje 



-j2-K{v'm+kTi ) /N 



where hi{n,m) — hi{nTs + [Lcp + m)T) is the sampled 
complex amplitude of the path, and k and v denote 
frequency and Doppler spread, respectively. Apparently, ICI 
is present due to the non-diagonal H/(ri). However, when 
fdTs < 0.1, the signal-to-interference ratio (SIR) is over 17.8 
dB [13], which enables us to discard non-diagonal elements 
of (n) with a negligible performance penalty. 

Por the comb-type pilot pattern, only pilot tones, denoted 
as yp{n) € C^^^, are extracted from yf{n). Then, approxi- 
mating H/(n) by a diagonal matrix, (|3]l is modified to 



yp{n) = Xp(n)hp(n) + np{n) 



(4) 



where Xj,(n) e C^^^ is a pre-known diagonal ma- 
trix, and hp{n) E C^^^ consists of the diagonal el- 
ements of H/(n). By denoting the instantaneous CPR 

as H{n,ra,k) = ^^"'^ /ii(n, m)e~-'-^^''"^'/^, we have 



\hp{n)]p = ;^ I]^j=o m, 6*^), where Op is the index 
of the p-th pilot tone. Hence, hp(n) is the time-averaging 
CPR during the n-th OPDM symbol. Besides, the noise term 
np[7i)^CN{0,allp). 

III. Maximum Doppler Spread Estimation 

Pirst of all, the least-squared (LS) channel estimation on 
pilot tones is carried out at the receiver, that is, 

hp,is(n) = Xp(n)"Vp('^) = hp(n) + Wp(n) (5) 

where h.p^is{n) E C^^^ is the LS-estimated time- 
averaging CPR, and Wp(rt) = ^p{n)~^np{n) is the noise 
term. Por equipowered and PSK-modulated pilot tones, 
Xp{n)"Xp{n) = Ip, therefore Wp(n) - CA/'(0, cr^Ip). 



By defining the Pourier transform matrix Fp E C 

[Fp,r]p,; = g-j^TrOpTi/N ^ jjjg 0-lag auto-correlation matrix is 

R,,(0) = E [hp,i,(n)hp,i,(n)^] = ^oFp,,DFf, + allp 

(6) 

where D — diag{(jf), I — 0, . . . , L — 1, and is 



^0 = ]^ E E JoC^^Um - q)T) (7) 

m=0 g=0 

Similarly, the /3-lag auto-correlation matrix of hp_;s(n) is 
R^(/3) = E [hp,is{n + I3)\ip,u{nf] = ^^Fp,,DF^, (8) 
where is 

= ]^ E E M^TTfdim - <z + /3(1 + r,p)N)T) (9) 

ni=0 

If the number of channel taps is known, the delay-subspace, 
denoted as Ut- s E C^^^, and the variance of noise af^ can 
be acquired with the eigenvalue decomposition. If unknown, 
however, the number of significant taps can be estimated by 
applying the minimum description length (MDL) algorithm 
[14]. Then, a Doppler dependent parameter is defined as 



^ / ||Ug R;,(/3)U.,,|||, _ 

^ \\lJ?MO)lJr.s-allL\\j. Co ^ ^ 

Through some manipulations (see Appendix [All, i^o and 
can be expanded into series and approximately expressed as 



eo = E^^--E 



oo 



k^O 

oo 



^/3 = E^fc^E 



fc=0 

1 



fe=0 



k\{k + l)l{2k + 1) 
k\{k + l)\{2k + 1) 



(11) 



2fe+2 



2ip 



2k+2-\ 



(12) 



where TV > 1, -0 = irfdNT, and = + r^p). 

Since the absolute values of Sk and tk decrease very fast, a 
finite number, say, K, of terms in ( ITTT i and (fTST i are sufficient 
to meet the accuracy requirement. So, with (fT0t-(fT2]i. we have 

El~o'(*fc-'?*^) =0 ^13) 

Let a; = — V'^ and Ck = 
Then, (T3[ is equivalent to 



2fc!(*:+l)!(2fe+l) 



V Ckx'' = 



(14) 



( fT4l l is a non-linear high-order polynomial equation. By 
resorting to Newton's method, its root, denoted as x*, is 
readily solved after numbers of iterations. Accordingly, fd is 
given by 



fd — 



ttNT 



(15) 



When the transceiver is moving, the path delays of the 
channel are slowly drifting [15] [8], which causes Fp;,- to vary 
and so does XJr.s- To accommodate the variation, the subspace 
tracking algorithm is adopted to automatically update Ut- s. 
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The proposed algorithm is summarized in the following. It is 
worth noting that its computation complexity depends on the 
QR-decomposition operation that is 0{P x L^). For sparse 
multipath channels, the complexity is moderate because L is 
quite small. Besides, the choice of K in (fT4] i can be made 
according to the tradeoff between the accuracy and complexity. 
After K is chosen, the convergence of Newton's method is 
quadratic along the number of iterations. In fact, numerical 
results show that less than 4 iterations are sufficient to achieve 
a precision of lO^"*. Moreover, for each iteration, the first 
differential is readily obtained thanks to the polynomial co- 
efficients, which reduces the complexity of Newton's method 
considerably. Finally, as is oscillatingly attenuating along 
/3, (fl4l i has multiple roots for large /3 0. In order to converge 
to the specific root, the initialization of Newton's iteration, 
therefore, should be carefully chosen according to (3. 

Initialize: (n = 0) 

Qo(0) = Q,5(0) = [lL^,Ol^,P-L„f 

Ao(0) = A0{O) = Op,L„, Co(0) = Ci3{0) = Il„ 
Run : {n — n + 1) 



Input : 


hp,is(n) 


Step 1 


: Updating eigenvalues of Rh(0): 


Zo(n) = 




Ao(n) = 


= aAo(n - l)Co(n - 1) + (1 - Q)Zo(n)Qo(7i - 1)" 


Ao(n) = 


= Qo(n)Ro(n) : QR- factorization 


Coin) = 


= Qo(n- l)^Qo(n) 


Step 2 


: Updating eigenvalues of Rh(/3): 


Z^(n) 






= aAp{n ~ l)C^(n - 1) + (1 - a)Zp{n)Qp{n - 1)' 




— Q^(n)R^(n) : QR- factorization 


Cp{n) 


= Q^(n-l)^Q^(n) 


Step 3 


: Estimating L and a'^: 



L = MDL (dzag(Ro(n))), = ^ Ep^L+i 1^(^)1 



Step 4 : Estimating r] according to (1101) : 

y ^;-j;'|[R0(n)],,,-*?.|^ 

step 5: Estimating fa by (fTIT) and (fT5l) . 
Remark: a is an exponential forgetting 
factor close to 1. Lm is the maximum rank to 
be tested. MDL{-) denotes the MDL detector. 
In simulations, we set a = 0.995, f3 — 1 , 
Lm ~ 10, K = 8, and the precision threshold 
of Newton iteration and the maximum number 
of iterations as 10""* and 4, respectively. 



IV. Simulation Results 

The performance of the proposed algorithm is evaluated for 
an OFDM system with BW = 12 MHz (T = 1/BW = 83.3 
ns), N = 1024, Lcp = 128 and P = 128. Two 3GPP E- 
UTRA channel models are adopted: Extended Vehicular A 
model (EVA) and Extended Typical Urban model (ETU) [16]. 
The delay profile of EVA is [0, 30, 150, 310, 370, 710, 1090, 

'For applied OFDM systems, I14t has only one negative root close to zero 
when l3 < 4. 



1730, 2510] ns, and its power profile is [0.0, -1.5, -1.4, -3.6, 
-0.6, -9.1, -7.0, -12.0, -16.9] dB. For ETU, they are [0, 
50, 120, 200, 230, 500, 1600, 2300, 5000] ns and [-1.0, -1.0, 
-1.0, 0.0, 0.0, 0.0, -3.0, -5.0, -7.0] dB, respectively. The 
classic Doppler spectrum, i.e., Jakes' spectrum [11], is applied 
to generate the Rayleigh fading channels. 

In Figlll the proposed algorithm is evaluated for fd = 200, 
400, and 600 Hz, respectively, within a wide range of SNR's. 
And two durations of observation (20 and 40 ms) are adopted 
to acquire the sample auto-correlation matrices. From the 
figure, it is evident that the performance of the proposed algo- 
rithm is robust when SNR > 5 dB, since the delay-subspace 
effectively reduces the noise disturbance. Furthermore, the 
estimation accuracy is rather high for fd > 400 Hz and the 
40 ms observation but somewhat deteriorated for 200 Hz and 
20 ms. It is accounted for the accumulation process of the 
sample matrices: the larger the maximum Doppler spread is, 
the faster the CFR updates; the longer the duration is, the more 
sufficient the sample matrices are. When fd is yet smaller, say, 
50 Hz, the proposed algorithm, nevertheless, may fail with 
the 40 ms observation due to the insufficient sample matrices. 
Besides, a simplified version of the proposed algorithm [17] 
still outperforms the CP-based algorithm [5]. 

The convergence of the proposed algorithm is verified with 
various Doppler's, i.e., fd = 600, 400, 200 Hz, when SNR = 
15 dB for EVA and ETU channels, respectively. As plotted in 
Fig|2] the estimated fd's converge after hundreds of samples. 
Furthermore, since the CFR updates faster when fd is larger, 
the convergence speed is faster for larger fd than smaller In 
addition, for all the curves drawn in Fig|2] the estimated values 
fluctuate around their true ones within a certain range, and the 
variations are larger for smaller fd's because of the sensitivity 
to the estimation error of rj in ( fTOl i when fdTg is small. If 
necessary, an averaging/smoothing window over the output of 
the proposed algorithm can be applied to supply a more stable 
estimation. 

V. Conclusions 

The maximum Doppler spread is crucial for adaptive strate- 
gies in OFDM systems. In this paper, we propose a subspace- 
based maximum Doppler spread estimation algorithm appli- 
cable to the comb-type pilot pattern. By tracking the drifting 
delay-subspace, the noise is greatly reduced. And by solving 
the polynomial equation with the Newton's method, high 
accuracy can be achieved with moderate complexity. The 
performance of the proposed algorithm is demonstrated to 
be robust over a wide range of SNR's and Doppler's by 
simulations. Besides, the proposed algorithm can be readily 
integrated into the channel estimators with the subspace tracker 
[8] [18], which lends a broad application promise to it. 

Appendix A 
Approximations of and 

The Maclakutin series of Jo{z) is [19] 
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Fig. 1. The performance evaluation for the proposed algorithms for EVA 
and ETU channels with different amounts of samples and Doppler's. 



Fig. 2. The convergence of the proposed algorithm when SNR - 
and l3 = 1 for EVA and ETU channels. 
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hence ^ is rewritten as 



m^ —0 mo—0 



{irfdT [(7711 - r7i2) + + np)N]f'' (17) 



Denoting tp = irfdNT and tp — /3(1 + Tcp), tk is further 
expanded as 

/ i2\h ^"-1 ^-1 2fe 



tk 



{kiNy 



E EE^2V=-x 



mi— m2— p—0 



0=0 
.2\k 



2k 



(fc!)2 



N-1 1 r , N-1 



N ^ ^ N ' N ^ ^ N ' 

7n-\ —0 



N ^ ^ N 

m2=0 



(18) 



When N ^ 1, we have 



, N-1 

_ V f— )« ~ ^ — V (~—)p-i ~ i 



1 ^-1 

1 1 ^ . 7712, 



i-iy- 



mi— ^ ' m2— 

Thus, tk can be approximated as 



q + l 



2k 



/c!(fc + 1)!(2A: + 1) 



i [(1 + <^)2'=+2 + (1 _ <^)2'=+2 _ 2(^2fc+2] (^9) 

Similarly, by letting /? = in ( fT9] l. (|2]l is expanded and 
approximated as 

^o = E^^«Efc!(fc + l),(2fc + i) ^20) 



fe=0 



fc=0 
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